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Abstract. This paper addresses some numerical and theoretical aspects of dual Schur domain 
decomposition methods for linear first-order transient partial differential equations. The spatially 
discrete system of equations resulting from a dual Schur domain decomposition method can be 
expressed as a system of differential algebraic equations (DAEs). In this work, we consider the 
trapezoidal family of schemes for integrating the ordinary differential equations (ODEs) for each 
subdomain and present four different coupling methods, corresponding to different algebraic con- 
straints, for enforcing kinematic continuity on the interface between the subdomains. Unlike the 
continuous formulation, the discretized formulation of the transient problem is unable to enforce si- 
multaneously the continuity of both the primary variable and its rate along the subdomain interface 
(except for the backward Euler method). 

Method 1 (d-continuity) is based on the conventional approach using continuity of the primary 
variable and we show that this method is unstable for a lot of commonly used time integrators 
including the mid-point rule. To alleviate this difficulty, we propose a new Method 2 (Modified 
d-continuity) and prove its stability for coupling all time integrators in the trapezoidal family 
(except the forward Euler). Method 3 (u-continuity) is based on enforcing the continuity of the 
time derivative of the primary variable. However, this constraint introduces a drift in the primary 
variable on the interface. We present Method 4 (Baumgarte stabilized) which uses Baumgarte 
stabilization to limit this drift and we derive bounds for the stabilization parameter to ensure sta- 
bility. Our stability analysis is based on the "energy" method, and one of the main contributions 
of this paper is the extension of the energy method (which was previously introduced in the con- 
text of numerical methods for ODEs) to assess the stability of numerical formulations for index-2 
differential-algebraic equations (DAEs). Finally, we present numerical examples to corroborate our 
theoretical predictions. 



1. INTRODUCTION 

Solving partial differential equations using domain decomposition techniques has been an active 
area of research for quite some time \22\ [26l I27j . Much of this research has addressed static 
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problems [9l [19], designing preconditioners [Ml [E], parallel implementation [9l [6], efficient 
computer implementation [9l[T2l[2T], solvers [25], etc. On the other hand, domain decomposition 
techniques for transient problems are not as mature as those associated with static problems. 
Some of the representative papers on domain decomposition methods for transient problems are 
[18\ [3 [6l [3 [T2| [5l [21]. A study of the choice of the interface boundary conditions has been done 
for second-order transient problems (i.e., structural dynamics) using the Newmark family of time 
integrators. For example, see Cardona and Geradin [3], Farhat et. al. [7], and Combescure and 
Gravouil [5]. All these papers deal with structural dynamics, which is a second order transient 
systems. 

In Reference |17j a coupling method has been proposed for first-order transient systems that can 
accommodate different time integrators and time steps in different subdomains. This was made 
possible by making the calculation of the interface Lagrange multiplier explicit. In this paper we 
restrict our study to time stepping schemes from the generalized trapezoidal family, and assume 
same time integrator and time step in all subdomains. This will allow the calculation of the interface 
Lagrange multiplier to be implicit, and thereby increasing the overall numerical stability. 

The main objective of this article is to present four variants of dual Schur domain decomposition 
method for linear first-order transient systems, and assess their stability. We restrict our study to 
conforming computational meshes. The stability analysis is based on the "energy" method |23[ll4j. 
The stability results for first-order systems (which are index-2 DAEs) presented in this paper are 
not a direct extension of results for ODEs or second-order systems (e.g., structural dynamics), and 
to the best of our knowledge, have not been reported in the literature. 

Remark 1.1. A thorough discussion of other computational aspects like preconditioners, parallel 
implementation issues (like scalability and speedup), direct and iterative solvers in the context of 
domain decomposition method is beyond the scope of this paper. These computational aspects have 
been addressed adequately in the literature and we have cited a few representative references to 
this end. These techniques can be generally employed without difficulty for the proposed domain 
decomposition methods. 

1.1. Main contributions. Some of the main contributions of this paper are 

• extension of the energy method for assessing stability to linear index-2 differential-algebraic 
equations, 
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• stability analysis of the d-continuity method, and demonstrated the instability of the d- 
continuity method under many common time integrators from the generalized trapezoidal 
family, 

• modified d-continuity method and its stability proof, 

• derived an upper bound for the user-defined parameter in the Baumgarte stabilization to 
ensure numerical stability, and also estimated the critical time step for the method, and 

• verified numerically that the optimal spatial convergence rate is unaffected by these domain 
decomposition methods. 

2. TIME CONTINUOUS GOVERNING EQUATIONS 

We will consider linear transient heat conduction as our model problem to illustrate various 
coupling methods. Consider a continuous domain with prescribed temperature on diQ and 
prescribed fluxes on ^2^^ (and for well-posedness we have diil. D = and difl U ^2^^ = dil.). 
The semi-discrete finite element equations for linear transient heat conduction can be written as 
(for example, see Hughes [13] ) 

(2.1) Mu + Ku = f, vte[o,r] 

where u{t) is the nodal temperature vector, t denotes the time, superposed dot denotes the time 
derivative, M is the capacity matrix, K is the conductivity matrix, and f{t) is the (prescribed) 
external heat source vector. The capacity matrix M is assumed to be symmetric and positive 
definite, and the conductivity matrix K is assumed to be symmetric and positive semidefinite. The 
following initial conditions and constraints complete the above differential system: 

(2.2) u{t = 0) = u'^^\ u\s,Q = Up{t) 

where is the prescribed nodal initial temperature, and Up is the prescribed nodal temperature 
on the boundary diQ. 

2.1. Decomposed problem. We decompose the computational domain into S subdomains fol- 
lowing a dual Schur formulation. This approach automatically implies that the equilibrium of 
the interface fluxes is enforced through the Lagrange multipliers. We assume that the kinematic 
constraints are linearly independent. We assume that the kinematic constraints are linearly inde- 
pendent and are compactly expressed using signed Boolean matrices. A signed Boolean matrix has 
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entries either —1, or +1 such that each row has at most one non-zero entry, see Reference jl7j . 
Note that by the construction of the signed Boolean matrix C^, its 1-norm is ||Cj||i < 1, which will 
be used in Section [S] Also note that if Cj / (which is the case in this paper) ||Ci||i = 1. 

Remark 2.1. Using signed Boolean matrices as defined above one can handle cross points (which 
are the points at which more than two subdomains intersect), and should be able to construct linearly 
independent kinematic constraints. For further details see Reference [T7] . 

The decomposition into subdomains of the differential system given by equation (j2.ip can be 
written as follows: Vt S [0, T] we have 

(2.3) M,u, + K,Ui = f, + C^^X Vi = !,•••, 5 

s 

(2.4) Y.C^Ui = 

i=l 

where A represents the vector of Lagrange multipliers, and Cj the signed Boolean connectivity 
matrices for perfect connection between compatible meshes. Note that the time continuous formu- 
lation is capable of enforcing the continuity of all the kinematic quantities and their time derivatives 
along the subdomain interface. However, this is not true of time discretized equations as we will 
see in the next section. 

Remark 2.2. For a derivation of equations (|2.3p - (|2.4p see Reference [n\ Appendix A]. 

It is important to note that the undecomposed problem (given by equation (j2.ip ) is a system of 
ODEs, whereas the decomposed problem consists of a system of ODEs given by equation ()2.3p . and 
a system of algebraic equations given by equation (j2.4p . Such a system of equations (that consists 
of a system of ODEs and a system of algebraic equations) belong to a broader class of equations 
called differential- algebraic equations (DAEs). For completeness, we present a very brief discussion 
of DAEs. 

An important special class is the semi-explicit DAE (also referred to an ODE with algebraic 
constraints), which can be mathematically written as 

(2.5) x = hi{t,x,y), = h2{t,x,y) 

It is easy to see that the governing equations of the decomposed problem, given by equations 
UM-^^M^ form a semi-exphcit DAE. 
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The index of a DAE is a measure to assess the difficulty to solve it numerically. Some of the 
common indices are - the differential index, local index, geometric index, tractability index, and 
perturbation index pQ. A popular index that is easy to compute is the differential index. For 
a semi-explicit DAE, the differential index can be defined as the minimum number of times the 
(algebraic) constraint must be differentiated to put the DAE in the standard ODE form 



using purely algebraic manipulations. The higher the differential index, the greater is the difficulty 
to solve the DAE numerically. For the DAE given by equations ()2.3p - (l2.4p . the differential index 
is 2. Note that the differential index of a system of ODEs is 0. For a thorough discussion of 
differential-algebraic equations refer to [H [16] . 

It is well-known that numerical solvers for ODEs may not work well for DAEs. In many cases 
there are stability and accuracy (e.g., drift in the kinematic constraints) issues. This fact has been 
discussed in a seminal paper by Petzold [20]. Time stepping schemes from the generalized trape- 
zoidal family have been developed primarily for numerically solving ODEs. The stability results for 
most numerical integrators that have been reported in the literature (for example. Reference |14j ) 
are appropriate for ODEs. The present paper addresses the stability of the time stepping schemes 
from the generalized trapezoidal family when applied to DAEs of the form given by equations 



The time interval of interest [0, T] is divided into N equal time steps of size At > and the 
equations are to be solved numerically at discrete instants of time tn = to + nAt with to = 
and ti\j = T. One of the popular ways to solve equation (|2.ip numerically is by employing a time 
stepping scheme from the generalized trapezoidal family. The numerical solution procedure for 





([13])- ([231). 



3. TIME DISCRETE EQUATIONS 



the generalized trapezoidal 



family can be written as: given (d^" '^\v^'^ obtain (d^"\i)(")) by 



solving 



(3.1) 



(3.2) 
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where < 7 < 1 is a user-defined parameter, and 

(3.3) « U{tn), V^^^ « U{tn), ^ f{tn) V?Z 



Some of the popular time integrators that belong to the generalized trapezoidal family are forward 
Euler (7 = 0), midpoint rule (7 = 1/2), and backward Euler (7 = 1). It is well-known that the time 
stepping schemes from the generalized trapezoidal family when applied to ODEs are unconditionally 
stable for 7 > 1/2 and conditionally stable for 7 < 1/2. For further details see Reference |14j . 

3.1. Decomposed problem. In the time discrete case of the decomposed problem two cases 
can be envisaged for enforcing the kinematic continuity. We can either prescribe continuity of 
temperature (which we refer to as <i-continuity) , or continuity of rate of temperature (which we 
refer to as i)-continuity) . These two kinematic continuity constraints can be mathematically written 
as follows. 

S 

(3.4) d-continuity: ^ Cidf ^ = Vn 

i=l 
S 

(3.5) iJ-continuity: ^ Cji^^-"'' = V?i 

i=l 

From a discrete point of view, we cannot, in general, enforce the continuity of both the temperature 
and the rates at the interface. (One exception, as we will show later, is the backward Euler 
with d-continuity in which we can satisfy both the discrete kinematic continuity constraints, see 
Remark [531) Based on the choice of kinematic continuity constraints we have two classes of domain 
decomposition method. We now present four different domain decomposition methods of which two 
methods employ d-continuity, the third method employs i)-continuity, and the fourth uses a linear 
combination of d-continuity and ?;-continuity. 



4. DUAL SCHUR DOMAIN DECOMPOSITION METHODS 

4.1. d-continuity domain decomposition method. This method can be written as: Vn = 
1, • • • ,N and Vi = I,-- - ,«S'; obtain ^d|"\ vj^^ A^"^^ by solving the following linear system of 
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algebraic equations. 

(4.1) M,^S") + KidS") = /S") + C,Ta(") 

(4.2) ) = d^-'^ + At ((1 - 7)^f + 7^S"^ 

(4.3) EC^<^"^=0 

4.2. Modified d-continuity domain decomposition method. This method is actually moti- 
vated by the stability analysis of the d-continuity domain decomposition method, which is presented 
in a later section. The modified d-continuity method can be written as: Vn = 0, • • • , — 1 and 
yi = 1, ■ ■ ■ , S; obtain ^d^^^^\ v^"'~^'^\ X^"''^'^^^ by solving the following linear system of algebraic 
equations. 

(4.4) M,^f + K,d^^^^ = + CTa(-+^) 

(4.5) df = df) + Ai^f+^); df = (1 - 7)df) + 7df+^) 

(4.6) X;C,d("+^)=0 

i=l 

If required, calculate v^"^^^ and A^"^^^ as follows (see Figured]). 

(4.7) VI = jv] + (1 - -yjv] 

(4.8) A("+^) = 7A("+^) + (1 - 7)A("+i+^) 

Remark 4.1. In the modified d-continuity domain decomposition method the equilibrium is enforced 
at time level {n + 7) and the kinematic constraints are enforced at time level {n + 1) (that is, at 
an integer time level). See Appendix for a numerical implementation procedure for the modified 
d-continuity method. 

Remark 4.2. One cannot employ the forward Euler (7 = 0) in the d-continuity and modified 
d-continuity domain decomposition methods. To see this, let us consider the d-continuity method. 

For the case 7 = 0, given d^" and f ^" , the calculation of d ."^ becomes explicit as we can 
directly evaluate the quantity using equation (|4.2p . The kinematic constraints given by equation 



3]) may not be consistent with the d-"^ computed from equation (j4.2p . We still need to find v^^'^ 

and A^"^ but we have only the subdomain equation (|4.1|) . which is also of the same size as the 
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modified tZ-continuity method 



cZ-continuity method 





H in 



tn+1 

X • 



[1 - 7)Ai 



-9 r„+2 



^n— 1+7 ^n+7 ^n+1+7 

Figure 1. This figure shows the main difference between the interpolations used 
in the d-continuity and modified d-continuity domain decomposition methods for 
calculating Lagrange multipliers and rates at various time levels. In the d-continuity 
method, one solves for the Lagrange multipliers (and rates) at the integral time levels 
tn- If one chooses, the Lagrange multipliers at the time levels tn+7 can be obtained 
as A*^'^"'"'''^ = (1 - 7)-^''"^ + 7A'^""''^\ which is illustrated in the middle figure. In 
the modified d-continuity method one solves for the Lagrange multipliers at the 
weighted time levels tn+7- If one chooses, the Lagrange multipliers at integral time 
levels can be obtained as A*^"^ = 7A'^"~"'^"'''^^ + (1 - 7)A'^""'"'^\ which is illustrated in 
the top figure. Similar explanation holds for interpolating the rates at various time 
levels. 

individual subdomain vectors vf^\ With forward Euler, the system given by equations (j4.ip - (j4.3p 
will be under- determined, and hence the problem is not well-posed. Similarly, one can show that 
the forward Euler is not compatible with the modified d-continuity method. 



4.3. u-continuity domain decomposition method. This method is based on the standard 
index reduction technique by analytical differentiation of the (kinematic) constraints. The method 



can be written as: Vn = I,-- - and Vi = 1, • • • , S"; obtain ,v^\x^"'^j by solving the 

following linear system of algebraic equations. 

(4.9) M,t;S") + ) = ) + C.^aW 

(4.10) ) = d^-'^ + At ((1 - 7)^f + 7^S"^) 

(4.11) i;c.^i"^=o 

Since the ?;-continmty method enforces the constraints on the rates, one may get significant drift 
in the original constraint (i.e., continuity of temperature along the subdomain interfaces) for larger 
time steps and for larger number of steps (i.e., longer simulation time). This irrecoverable drift 
in the original constraint makes the ir-continuity method undesirable in such situations. One of 
the popular ways to control drift in the constraints is Baumgarte stabilization, which has been 
proposed in the context of multibody dynamics [2]. With a slight modification, one can extend 
Baumgarte stabilization to index-2 DAEs first-order transient systems (which is the case in this 
paper). Following this approach, we construct a new domain decomposition method, which is 
outlined in the next subsection. 

4.4. Baumgarte stabilized domain decomposition method. This method enforces a linear 
combination of the d-continuity and i>-continuity as the kinematic constraints. The method can 
be written as: Vn = 1, • • • , N and Vi = 1, • • • , S*; obtain (^d["'\ vf^\ A^"^^ by solving the following 
linear system of algebraic equations: 

(4.12) M,t;f ) + Kid^^ = + C^aW 

(4.13) dt^ = df^ + At ((1 - ^)vt'^ + 7^"^) 

s s 

(4.14) E'^^-!"^ + |7E^^^^^ = 

1=1 1=1 

where a > is a dimensionless user-defined parameter, which we will refer to as the Baumgarte 
parameter. As we will see in Section O the choice of a will effect the accuracy (the drift in the 
original constraint, that is, continuity of temperature along the subdomain interface) and stability 
(i.e., the critical time step). Typically, the larger the Baumgarte parameter, the smaller will be the 
drift in the original constraint, but also smaller may be the critical time step (which may depend 
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on the choice of 7) . In Section [5] we derive bounds on the choice of a that ensures stabihty of the 
domain decomposition method. 

Also, note that, unhke the d-continuity and modified d-continuity methods, one can employ the 
forward Euler (7 = 0) under the Baumgarte stabilized domain decomposition method. 

In the remainder of the paper we will address the following questions by providing rigorous 
mathematical justifications. 

• Under what conditions do the d-continuity and t)-continuity domain decomposition methods 
give stable results? 

• Do the algebraic equations alter the stability of the numerical time stepping schemes from 
the generalized trapezoidal family? If so, which of these time stepping schemes are stable? 

• Do any of the time-stepping schemes from the generalized trapezoidal family achieve con- 
tinuity of both temperature and rate of temperature along the subdomain interface? 

Remark 4.3. Before we perform stability analyses of the methods presented in the previous sec- 
tion, we briefly comment on the implementation of dual Schur domain decomposition methods. 
Various ways of implementating dual Schur domain decomposition methods have been proposed in 
the literature. A thorough description is beyond the scope of this paper. But, herein, we just cite 
few representative papers that can be consulted for effective implementation of dual Schur domain 
decomposition methods. 

The FETI method is a popular way of implementing the dual Schur domain decomposition method 
which has been shown to possess good convergence and scability properties [9l [H [8] . An interest- 
ing approach (in the context of structural dynamics) has been developed in \12\ I21j . Any of the 
aforementioned implementation methodologies can be employed for the methods presented in this 
paper. 

5. STABILITY ANALYSIS 

In this section we assess the stability of the d-continuity, modified d-continuity, ^-continuity and 
Baumgarte stabilized domain decomposition methods using the "energy" method [23\ I14j . Note 
that the energy method provides sufficient conditions for stability. However, in many instances the 
obtained bounds are quite sharp. For example, in the case of ODEs, the obtained stability condition 
using the energy method (that is, the critical time step) has turned out to be both necessary and 
sufficient |14j . 
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Let M"* denote the standard m-dimensional Euclidean space. We say the vectors aj^"-* G M"* (n 
0, 1, 2, • • • ) are bounded Vn if there exists a positive constant C independent of n such that 



(5.1) 



|^;(")||<C7 Vn 



where || • || is some convenient norm defined on (for example, say the 2-norm). Note that in 
finite dimensional vector spaces all norms are equivalent |13j . 

We employ the following notation for the jump and average operators over a time step: 

1 



(5.2) 



X 



(n) 



^{n+l) _ ^(n) 



/ 2 



X^n+l) _|_ ^(n) 



It is easy to show the following identities: 
(5.3) (1 - 7)a;(") + 7a;("+i) 

and for any symmetric matrix S we have 
(5.4) 

For convenience we define the matrix Ai as 
(5.5) 



X 



in) 



+ 





1 






~ 2 





1 



Ai := M, + (^7 - -j AtK, 

We choose the time step Ai in such a way that it satisfies the stability requirements for all 
individual unconstrained subdomains (i.e., A = 0). That is. 



(5.6) 



At < min(At5™,-- - ,Atf') 



and the critical time step for subdomain i is given by 

2 



(5.7) 



At: 



crit 



ujf^--{l-2'y) < 7 < 1/2 

+00 1/2<7<1 



where u;™^^ is the maximum eigenvalue of the generalized eigenvalue problem 
(5.8) iOiM,cl)i = Kicjji 

Note that all the eigenvalues tOi are real and non-negative [H]. For the chosen time step as described 
above, the matrices Aj (i = 1, • • • , S*) (defined in equation (j5.5|) ) are positive definite. Also, note 
that the matrices Ai are symmetric. 
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5.1. Stability of the d-continuity method. Equation (j4.ip implies 

(5.9) MiV^/"^^^ + = + C;^A("+^) 

where A("+'^) := (1 - 7)A(") + 7A("+i); /["+^^ := (1 - 7)/S"^ + 7/S"+^^ and similar expressions for 
and 



(5.10) 



Using equations (j4.2p . (|5.3p and (j5.5p : the above equation can be written as 



At 



d 



in) 



+Ki{4")} = /;"+^)+cTA("+^) 



Premultiplying the above equation by the vector |<i^"^|, using the identity given in equation (j5.4p . 
and then summing over all S sub domains; we get 



2At 



j=l i=l 

For stability analysis, one assumes the externally applied forces to be zero (i.e., /^"^ = 0; Vi 
1, • • • ,5; Vn) [H]. Thus, we have 



s 



i=l 



2At 



+ 



E{4"'}V,{4">} = E{«!"'}'c?A 



(n+7) 



i=l 



j=l 



By invoking the fact that the matrices Ki (i = 1, • • • ,5*) are positive semidefinite, we conclude 
that 



1=1 



2At 



A,-d^) 



< 



j=i 



i=l 



Using the linearity of the average operator (which allows us to interchange the summation and 
average operation), and the d-continuity given by equation (j4.3p : we conclude that 

^ 1 " 
^ 2At 



(5.14) 



i=l 



d,(") Aid,(") 

7. l' 1. 



< V?i 



Using the definition and linearity of the jump operator, we conclude 



s 



(5.15) 



i=l i=l 



(0)^ A.JO) 

i 



i=l 



Since Vi = 1, • • • , S the initial vectors df'^ are bounded and the matrices Aj are positive definite; 
we conclude that the vectors d^"^ [i = 1, • • • ,5) are bounded Vn. This implies that, from the 



trapezoidal equation (j4.2p . the vectors ■u^""'''^^ := (1 — 'y)vf^' +71'^"^^'' {i = 1, ■ • • ,S) are also 
bounded Vn. 
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We now show that the Lagrange muhiphers A*-""'"'''-' are bounded Vn. There are several ways in 
proving this resuh. Herein, we show using the Schur complement operator (which is widely used 
in computer implementations, see Appendix). An alternate derivation is presented in subsection 
15.31 for proving a similar result. Under zero external force the subdomain governing equation (j4.ip 
implies 

(5.16) M,t;;"+^)+K,dS"+^) =CfA("+^), Vi 

Using 1^;"+^) = (df - dS")) /Ai (see equation gj)) and d^^^^ = (1 - j)df'^ + the 
above equation can be written as 

(5.17) = {Mi - (1 - j)AtKi) df ) + AtM7^C^A("+^) 

(5.18) where Mi := Mi + jAtKi 

By premultiplying equation (I5.17P with Cj, then summing over the number of subdomains (i.e., the 
index i = 1, • • • ,S), and using the kinematic constraint for d-continuity method (equation (j4.3p ): 
we get 

(5.19) A("+^) = -^G-i ' - (1 - 7)AtK.) dS") 

i=l 

S 

(5.20) where G:=^C,M7^Cf 

1=1 

Since the vectors d^"^ are bounded Vn, from equation (|5.19p . we conclude that the Lagrange mul- 
tipliers A^""'"'^^ are also bounded Vn. 

Remark 5.1. Since the matrix Mi is positive definite (and hence invertible), all the steps in 
equation (jS.lTp are valid. Also, since the constraints are assumed to be linearly independent, one 
can easily show that the matrix G (which is sometimes called the Schur complement operator) is 
positive definite (and hence invertible). Therefore, all the steps in equation (j5.19p are valid. Also 
see Appendix. 

But the proof of stability is not yet complete as we have not said anything about the boundedness 
of v^"'^ and A^"^ and these quantities are used to advance the solution. Before we comment on the 
boundedness of rates and Lagrange multipliers at integer time levels t„, we state and prove the 
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following general result, which will be useful in assessing stability of coupling algorithms. We employ 
the standard notation used in mathematical analysis. To avoid ambiguity, let N := {0, 1, 2, • • ■ }. 

Proposition 5.2. Let be a sequence of real numbers and let := (l—'y)s^^'^+js^'^^^\ 

whereO < 7 < 1. If the element s^^^ andthe sequence {s^'^^"'^) are bounded, then for any j > 1/2 
the (original) sequence {s^"'^)^^^ is also bounded. 

Proof. We split the proof into two different cases: 7=1 and 1/2 < 7 < 1. 

First case (7 = 1): The proof is trivial as in this case s^^'^'^^ = s^"'~^^\ The sequence {s^^^)^^^ 
can be obtained by augmenting the element s^^^ at the start of the sequence {s^^^"'^) Since 
the element s^^^ and the sequence are bounded, we conclude that the original sequence 

(s^"'^) is also bounded. 

Second case (1/2 < 7 < 1): Construct a new sequence 

(c^''^)„gN s^c^ that c(") = If 
('^''"^)nGN bounded then {c^"'^)^^^ is also bounded and vice-versa. We now prove the proposition 
by the method of contradiction. 

Assume that the sequence (c^"''')^gj^ is unbounded. Then we can find a strictly increasing un- 
bounded subsequence (nk € N, < ni < n2 < • • • ) such that 



(5.21) 
(5.22) 



< c("i) < c("2) < ,00 and 



Since the sequence i^^"''^'^^) bounded, we can find < Af E R such that 



(5.23) 



< M Vn 



Since the subsequence (c*-"''°^)^gp^ is strictly increasing and unbounded, we can find an element of 
this subsequence such that 

1 



(5.24) 



-M, p G N 



27 - 1 

Since 1/2 < 7 < 1 (and, therefore 7 > 27 — 1), from equations (j5.23p and (|5.24p we have 



(5.25) 



7c 



(rip) 



> (27 - l)c("f) > M > 



Jn+y) 



> Vn 



Using the definition of s^'^p ^^"'^ we have 



(5.26) 



g{"p-l) 


1 




" 1-7 
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By using the triangle inequality (and noting that c^"'''^ := Is^"*"^!) we conclude that 



(5.27) 



.(rip-l) 



> 



1 



1-7 

By using equation (j5.25p we obtain 



{np-l+7) 



1-7 



7C 



irip) 



(5.28) 



> 



1-7 



^c("p) - (27 - l)c("'') 



which is a contradiction as it violates equation (j5.22|) . This implies that the sequence (c^"^)^gj^ is 



bounded, and so should be the sequence (s 



in)] 



'neN' 



This completes the proof. 



□ 



Remark 5.3. The result proved in Proposition 1 5. S\ in general, cannot he extended to < ^ < 1/2. 
Counterexamples for the cases < 7 < 1/2 and 7 = 1/2 are given in figureslM respectively. 
As discussed in Remark the forward Euler (7 = 0) is not well-posed under the d- continuity 
method. This implies that, under the d-continuity domain decomposition method, many of the 
popular time stepping schemes from the generalized trapezoidal family (e.g., the forward Euler and 
midpoint rule) are unstable. The midpoint rule (7 = 1/2) is unconditionally stable for linear first- 
order ODEs. But for linear index-2 DAEs, the midpoint rule can he unstable (and is right on the 
boundary of the instability region). 

Returning to the proof of stability of the d-continuity domain decomposition method, we have 
proved that v^^^'^'^ and A'-""'"^-' are bounded Vi and Vn. Applying Proposition 15.21 for individual 
components of the vectors i?-""'''^^ and A^'^''''^^ we can conclude that for 1/2 < 7 < 1 the vectors 
and A^"^ are bounded Vi and Vn. As shown mathematically in Remark 15.31 the quantities i^^"^ 
and A^"^ may not be bounded when 0<7<l/2. In a later subsection we will present physical 
systems that, in fact, exhibit this kind of (numerical) unbounded behavior for 7 < 1/2 under the 
d-continuity method. This completes the stability analysis of the d-continuity method. 

5.2. Stability of modified d-continuity method. The majority of the proof for this method 
is identical to the initial part of the proof for the d-continuity method (presented in subsection 
15. ip up to the step where we deduced that the vectors A^""'"'^^ and v^[^^'^'^ are bounded. The only 
additional thing we have to prove is the boundedness of the vectors A'^"'-' and t;^"^ , which is a direct 
consequence of the triangle inequality. To wit, using equation (|4.8p (and also see Figure [1]) we have 



(5.29) 



A(-) 



^Xin-l+l) + (1 _ ^)A("+T) 
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<7 



^(ri-1+7) 



+ (1-7) 



A("+'y) 



7 < 1/2 7 7 7 




> n 



Figure 2. A counterexample for the case < 7 < 1/2. The figure presents an 
example in which the (constant) sequence [s^"^^^') is bounded but (^s^^'^yi,^^ is 
unbounded. The elements of {s^'^~^'~^^)n&n are indicated using filled squares, and 
those of (s('"))„gN are indicated using circles. Recall that s^'"^'''-* := (1 — 7)s('") + 
7s("+i), Vn. 



n 





1 


2 


3 


4 ... 


s("+l/2) 


+1 


-1 


+1 


-1 


+1 ••• 


sin) 





+2 


-4 


+6 


-8 ••• 



Figure 3. A counterexample for 7 = 1/2. The sequence {s^^^^^'^^) > s("+i/2) — 
i^s^n) _|_ g{n+i)^ ig bounded but the sequence (•s^"^)„gf^ is unbounded. 

Since the vectors A^""*"^^ are bounded V??,, from the above equation it is evident that the Lagrange 
multipliers at integral time levels A^*^-* are bounded Vn. Using a similar reasoning, since the vectors 
^("+7) g^j^g bounded Vn, the rates at the integer time levels are also bounded. (Note that the 
vectors d^^^ are also bounded, which has been proven in subsection 15.11 ) This means that under 
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this method all quantities of interest are bounded. This completes the stability analysis of the 
modified <i-continuity method. 

Note that, as discussed in Remark 14.21 the forward Euler (7 = 0) is not well-posed under the 
modified d-continuity method. All other time integrators from the generalized trapezoidal family 
(0 < 7 < 1) can be employed and are stable under the modified d-continuity domain decomposition 
method. 



Remark 5.4. For the d-continuity and modified d-continuity domain decomposition methods, with 
backward Euler (7 = 1) one can achieve the continuity of both temperatures and temperature rates 
along the subdomain interface. To wit, the d-continuity constraints (which are basically the conti- 
nuity of temperatures along the interface, and are given by equation (j4.3p or (j4.6p ) imply 

s 

(5.30) J^C, 



i=l 



d 



(n-l) 



0, Vn 



For backward Euler we have 
subdomain interface. 

(5.31) 



(n-l) 



(n) 

Atv- . This implies the continuity of rates along the 



i=l 



5.3. Stability of the t>-continuity method. We follow a similar procedure employed in sub- 
section 15.11 Using equations ()4.ip - (|4.2p (and as usual neglecting the external forcing function for 
stability analysis [H]; i.e., f^^^ = 0) one obtains 



(5.32) Yl 



1 



+ 



i=l 



i=l 



Since the matrices Ki {i = 1, - ■ ■ ,S) are positive semidefinite, and At > 0; we conclude 



(5.33) 



^ 2 

j=i 



i=l 



Using the D-continuity given by equation (j4.1ip we conclude that 

s 

(5.34) 



^ 1 



i=l 



< Vn 



Using the definition and linearity of the jump operator, we conclude 

s 



(5.35) 



s s 



i=l 
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Since Vi = 1, • • • , S the initial vectors vf'^ are bounded and the matrices Aj are positive definite; 
we conclude that v["'^ are bounded Vn. This implies that, from the trapezoidal equation (|4.10|) . 



(n) 



is bounded. 



One can further show that, under the i)-continuity method, the jump in the Lagrange multipliers 
A^"^ is also bounded. There are several ways to prove this result, which was the case even 
under the d-continuity method for proving the boundedness of the Lagrange multipliers A^"'"'''^^ 
(see subsection 15. ip . As mentioned earlier, herein we take a slightly different approach (than the 
one in subsection 15. ip . and which is also applicable for the Baumgarte stabilized method for showing 
a similar result. To this end, we start with 



(5.36) 



Mi 



d 



in) 



c 



Vi 



(Note that in the above equation we have used the fact that the external force on all subdomains 
is zero, /^"^ = 0.) By premultiplying the above equation with CiM^ ^ , and then summing over 
the number of subdomains, we get 



(5.37) 



Mi 



V 



in) 



+ Ki 



d 



in) 



where the matrices Mi and G are defined in equations (|5.18p and (|5.20p . respectively. Since the 



vectors 



d 



in) 



and i» -"^ (and hence 



,("■) 



) are bounded Vn, from equation (j5.37p . we conclude that 



A(") 



is bounded Vn. 



Xin) 



is bounded), in the step just 



Remark 5.5. In order to show the above result (that the vector 
above equation (j5.37p . one can premultiply with any matrix of the form CiHi where Hi is some 
positive definite matrix. Since Mi ^ is a positive definite matrix (as discussed in Remark \ 5. 1\ that 
the matrix Mi is positive definite), the choice Hi = Mi ^ is valid. We have chosen this particular 
choice as to be able to use some of the earlier results (e.g., the matrix G is invertible), and to avoid 
introducing additional notation. 



For the i7-continuity domain decomposition method, solely based on the energy method, one 
cannot infer the boundedness of the quantities d-"^ and X^'^\ Also, there can be drift in the 
continuity of temperatures along the subdomain interface. That is, X^^i Cidf^^ ^ 0. This drift 
may grow over time due to round-off errors. As discussed earlier, one of the ways to control the 
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drift is to employ the Baumgarte stabilization. In the next subsection, we discuss the stability of 
the Baumgarte constraint stabilization in the context of index-2 linear first-order transient systems. 



Remark 5.6. Note that, unlike the d-continuity and modified d- continuity methods, one can em- 
ploy the forward Euler method (7 = 0) in individual subdomains under the v-continuity domain 
decomposition method. 

5.4. Stability of Baumgarte stabilized domain decomposition method. We start the sta- 
bility analysis by rewriting the kinematic constraints. To this end, equation (j4.14p implies that 



(5.38) 



S 

i=l 



An) 



+ 



a 
At 



(n) 







By using the kinematic constraints given by equation (j4.14p and identity ()5.3p we get 

s 



(5.39) 



i=l 



a 



+ a 



where the (dimensionless) parameter a* is introduced for convenience, and is defined as 



(5.40) 



a* :=l + a[j-- 



We now rewrite the subdomain equation in a manner similar to what we did above for the 
kinematic constraints. Equation (j4.12p implies that 



(5.41) 



d 



(n) 



Xin) 



By using the trapezoidal equation ()4.13p and employing the identity (j5.3p we get 



(5.42) 



Mi 



in) 



+ AtK, 7 



in) 



Premultiplying both sides of the above equation by ( a* 
over all of the subdomains we get 

s 



An) 



+ a 



1^^"^}) then summing 



E 



i=l 



a 



An) 



+ a 



in)\ 



(5.43) 



E 

i=l 



a 



V 



(n) 



-|- a 



V 



(n) 



+ AtKi (^(7 - 1/2) 



V 



a 



(n) 



+ 



in) 



i=l 



»!"'})) 
{4"'}) 



+ a 
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For the Baumgarte stabilized domain decomposition method, using equation (j5.39p . we conclude 
that 

s 



(5.44) 



(n) 



+ a 



i=l 



'"I \ 



Mi 



An) 



+ MKi (7 - 1/2) 



V 



(n) 



+ 







Using the definition of a* (equation (|5.40p ). the above equation can be rewritten as 



[vT^] +a{^(")}) V. [vt^] +AtY: [vl'-^K, ((7-1/2) 



1=1 



i=l 



(5.45) +aAtj2 ((7 - 1/2) [^S"^] + ((7 - 1/2) 

i=l 

Since a > 0, At > and l^j is positive semidefinite, we conclude 



{^-}) 



(5.46) 



{n) 



+ a {vt^]f [vf^] +Atj2 [^S"^] ^ ((7 - 1/2) [^f + {vf }) < 



i=l 



By invoking the symmetry of Mi, using equation (j5.40p . and rearranging the terms we get 

S rr, S 



(5.47) (aM, + AtK,0{^l"^} + EH"^ 



< 



i=l 



where the matrix A,- is defined as 



(5.48) 



Ai := a* Mi + At 7 



If the matrices Ai {i = 1, ■ ■ ■ , S) are positive semidefinite (later we will obtain sufficient condi- 
tions for the matrix Ai to be positive semidefinite) then equation (15.470 implies that 

(5.49) EK"T iaMi + AtKi)[v[''^]<0 

1=1 



By invoking identity (j5.4p we have 

s 



(5.50) 

This implies that 

s 

(5.51) Y 



E 

i=l 



) " {aM, + AtKi) ^ 



(n) 



< Vn 



(n+l) 



5 



,(0) 



i=l «=1 

Since a > 0, Mi is positive definite, and Ki is positive semidefinite; the matrix aMi + AtKi is 

positive definite. This implies that the vectors ■u^"^ are bounded Vn as the (initial) vectors vf'^ 
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are bounded. From the trapezoidal equation (|4.13p . one can conclude that the vectors 



are 



bounded Vn. Using the same derivation as presented under the t;-continuity method, one can 



conclude that the vector 



is bounded Vn. 

Furthermore, one can easily show that under the Baumgarte stabilized method the drift in the 
original constraint (i.e., continuity of temperature along the subdomain interface) is also bounded 
(but may not be zero), which is not the case with the v-continuity method. To show this (for 
convenience) let us work with 1-norm, and the result (boundedness) will hold for other norms also 
(because of equivalence of norms in finite dimensional vector spaces). We start with the 1-norm of 
the drift in the original constraint, and use equation (|4.14p and the triangle inequality: 



(5.52) nj^c^dr'n^^iitcr'ns^h^A 



a — ' a 

i=l i=l 1=1 



Now using the definition of matrix-norm, and the boundedness of v^""^; equation (j5.52p can be 
written as 



S3 3 

(5.53) llEC.df^ll < ^^IIQIIIKII < < C Vn 

i=l i=l i=l 

where C > is some constant independent of n, and we have used the fact that the 1-norm 
He'll!! — 1- (Also recall that if d ^ 0, which is the case in this paper, ||Ci||i = 1.) Equation 
(j5.53p implies that the drift in the original constraint is bounded Vn. 

Similar to the i^-continuity method, even under the Baumgarte stabilized method, one cannot 
infer anything about the boundedness of d^"'^ or Lagrange multipliers A^"^ solely based on the 
energy method. But, for many practical problems one often gets bounded numerical results for 
these quantities. The boundedness of the drift in the original constraint is the only additional (but 
very important) feature that the Baumgarte stabilized method has compared to the D-continuity 
method. 

We will now obtain sufficient conditions for the matrix Aj to be positive semidefinite (which we 
have assumed in the stability analysis of this method, see the line below equation (j5.48p ). Using 
the eigenvectors of the generalized eigenvalue problem ()5.8p as the basis, the matrix Ai can be 
diagonalized and obeys the similarity 

(5.54) Ai ~ a* I + At (^7 - fti 
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In the above equation, the symbol "~" denotes similarity of matrices, and the matrix fij is defined 

as 

(5.55) rii = diag {{ui)^ , • • • , {uJi)J 

where n is the size of the (square) matrix Af, and {uji)^ ■ ■ , (wi)„ are the eigenvalues of the gener- 
alized eigenvalue problem (15. Sh . (It is well-known that similar matrices have the same characteristic 
polynomial, and hence the same eigenvalues [13j.) 

Since Mi is positive definite and Ki is positive semidefinite, the eigenvalues (coi)-^ , • • • , (i-fj)„ are 
all non-negative. This implies that (using equations (|5.40p and (j5.54p . and noting that At > 0) the 
matrix Aj is positive semidefinite unconditionally for 7 > 1/2. For 7 < 1/2, sufficient conditions 
for the matrix Aj to be positive semidefinite are 

1 



(5.56) a* > =^ a< 



I7-I/2I 



(5.57) At < Atf = ^ = ^ — 

^ ^ - ' I7 - l/2|u;™'^^ (l-27)u;fi^^ wf^^^ 

Remark 5.7. Note that in equation ([5371) Atf'^ >0 as a* >0 and w™^'' > 0. 

These results imply that the use of Baumgarte stabilization decreases the critical step for subdo- 
main i by a/ivf^^^ relative to the unconstrained case. Summarizing the results, sufficient conditions 
for the stability of Baumgarte stabilized domain decomposition method are 



(5.58) < 7 < 1/2 



a < 



YPTJ2\ 



2 



(5.59) 1/2<7<1 < 



a < +00 
Atf^ = +00 

For a quick reference, the stability results derived in this section are summarized in Table [TJ In 
the next section we verify the theoretical predictions using numerical experiments. 

6. NUMERICAL RESULTS 

6.1. Numerical verification of stability results: Split-degree of freedom. A simple problem 

that is commonly used to study the performance of coupling algorithms is the split-degree of freedom 

system, which can be obtained by splitting a single degree of freedom into two single degree of 

freedom problems but subjected to kinematic and dynamic constraints, see Figure [H 
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Table 1 . A summary of stability results 



Method 



Stability characteristics 



d-continuity Stable for 1/2 < 7 < 1. 

Modified d-continuity Stable for < 7 < 1. 
i>-continuity Drift on the interface. 

Baumgarte stabilized Unconditionally stable for 1/2 < 7 < 1. Stability 

conditions for < 7 < 1/2 are given in Eq. (j5.58p . 



kA 



+A -A 



^ 7 7 ^ 7 



UB 



TUB 



Figure 4. Problem definition: Split-degree of freedom system 



The time discrete governing equations for the subdomains A and B consists of 



(6.1) 



ruAv'^^^ + fc^dj) = +A(-), ruBvP + ksdP = -A^ 



and along with the d- or i)-continuity (kinematic) constraints, which are respectively given by 



(6.2) 



- = 0, t;? ) - = 



In addition, we have the equation for the generalized trapezoidal family: 



(6.3) 



<B = <B'^ + At((l-7)^i,V+7<,B 



(n-1) , (n) 



In this subsection, using the split-degree of freedom, we numerically verify the derived theoretical 
stability results for the d-continuity and modified d-continuity methods. The system properties are 
taken as ttia = ms = 1, kA = W and fee = 1. The initial temperature is taken to be unity. The 
time step is taken as At = 0.01 s. To show all the aspects of the stability results, we consider two 
different choices of the trapezoidal parameter 7. 

In the first case we chose 7 = 1/4. For this choice, the critical time steps for (unconstrained) 
subdomains A and B are 0.4 s and 4 s, respectively. The obtained numerical results using the 
d-continuity method for the rate of temperature, Lagrange multiplier, and temperature are shown 
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in Figures [6l [7] and [8l respectively. As one can see from the figures, the obtained numerical results 
are unstable in accord with the theory presented in subsection 15.11 

For the same case (7 = 1/4), in Figure Owe have shown the rate of temperature and Lagrange 
multiplier that are obtained using the modified d-continuity method. Both the rate of temperature 
and Lagrange multiplier are all bounded even at the integer time levels, which agrees with the 
theory presented in subsection 15. 2i 

In the second case we chose 7 = 3/4. This time stepping scheme is unconditionally stable 
when applied to ODEs. Even when applied to DAEs, as shown in earlier in this paper, the time 
stepping scheme should be unconditionally stable with respect to all quantities. The obtained 
numerical results for the rate of temperature and Lagrange multiplier are shown in Figure [T0\ and 
the temperature is shown in Figure [TTl Even for this case, the obtained numerical results are stable, 
which agrees with the theoretical predictions. 

6.2. One- and two-dimensional problems, and multiple subdomains. The model problem 
is transient linear heat conduction, and the governing equation can be written as 



where u{x,t) is the temperature, V denotes the spatial gradient operator, f{x,t) is the volumetric 
source, x is the position vector (in 2D, x = {x,y)'^), k is the coefficient of conductivity, p is the 
density, and Cp is the coefficient of heat capacity. Li all our numerical simulations we have employed 
the standard Galerkin finite element formulation for the spatial discretization |14j . 

We consider two test problems, one each in ID and 2D. For the ID problem, a domain of length 
L = 2.0 m is divided into two equal subdomains, and each sub domain is divided into 10 equal 
linear finite elements. For the 2D problem, a square of size 2 m x 2 m is divided into four equal 
subdomains, and each subdomain is divided into 10 x 10 square elements (see Figure [5]). 

For both the test problems we assume homogeneous Neumann boundary conditions, and zero vol- 
umetric source (i.e., / = 0). The analytical solutions for ID and 2D test problems are, respectively, 
given by 



(6.4) 



pcpii — V ■ (kVu) = f 




■2 1 



(6.5) 



(6.6) 
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The problems are driven by prescribed non-zero initial conditions, which are obtained by evaluating 
the analytical solution (given in equations ()6.5p - (j6.6p ) at time t = 0. We use the standard linear 
and bilinear finite elements for ID and 2D test problems, respectively. We assume the material 
properties to be /c = pcp = 1. 

In Figure [12] we have shown the rate of temperature from the d-continuity method using 7 = 1/4. 
As predicted by the theory, even for the ID test problem, the numerical results are unstable. In 
Figure [13] we have shown the temperature and rate of temperature under the d-continuity using 
7 = 3/4, and as predicted by the theory the numerical results are stable and match well with the 
analytical solution. In Figure [T^ we have plotted the temperature obtained using the modified d- 
continuity method for both 7 = 1/4 and 7 = 3/4. For both the values of the Baumgarte parameter, 
the numerical results are stable and match well the analytical solution. 

In Figures [15] and [16] we have shown the numerical results for the 2D test problem from the 
d-continuity and modified d-continuity methods, respectively, for 7 = 3/4 and time step At = 
0.001 s. For both these methods, as expected, the obtained numerical results matched well with 
the analytical solution. 

6.2.1. Spatial numerical h-convergence analysis. From the standard approximation theory [3], it 

is well-known that the theoretical rate of /i-convergence for the above test problems (which have 

smooth solutions) using linear finite elements is 2. For a robust converging domain decomposition 

method, one would prefer to have the numerical rate of convergence to be the same as for the 

undecomposed problem, which for the chosen test problems is 2. The numerical spatial convergence 

analysis for linear transient heat conduction is typically done either by simultaneously decreasing 

the time step proportional to /i^ or by choosing a very small time step so that the error due to 

the time discretization is negligible. In our numerical simulations, we have employed the later 

approach, and have taken the time step to be At = 0.00001 s. We now investigate the numerical 

rate of /i-convergence for the d-continuity and modified d-continuity methods. 

The numerical convergence analysis is performed at time level t = 0.01 s using a hierarchy of 

uniform meshes. A typical 2D mesh along with subdomains is shown in Figure [5] In Figure \T7\ 

the numerical convergence results for ID and 2D test problems for the temperature and rate of 

temperature using the d-continuity method with the trapezoidal parameter 7 = 3/4 are shown. 

In Figure [T8l we have shown the numerical convergence results for the two test problems for the 

25 



temperature under the modified d-continuity method using 7 = 1/4 and 7 = 3/4. Under both 
the d-continuity and modified d-continuity methods and for both the test problems, the obtained 
numerical rate of spatial convergence is 2, which implies the methods are convergent schemes as 
predicted by the theory. 

6.3. Baumgarte stabilized domain decomposition method. In this subsection we study the 
numerical performance of the Baumgarte stabilized domain decomposition method, and compare 
numerical results with the theoretical results derived in subsection 15. 4i 

Consider linear transient heat conduction in a one-dimensional bar of length L = 2 m with 
material properties to be = 1 and pcp = 1, and the initial temperature to be unity. At the left 
end of the domain the fiux is zero, and on the right end we prescribe a constant temperature of 
zero. The analytical solution is given by 



4^ (-1)" r {2n + lY^:H^ \{2n + l)^x 



■ exp 

I zn -I- II 

n=0 



(2n + ly-nH 



cos 



2L 



4L 

The computational domain is divided into two equal subdomains, and each subdomain is divided 
into 10 equal finite elements. The left subdomain (which has Neumann boundary condition on 
its left end) is denoted as subdomain 1, and the right subdomain (which has Dirichlet boundary 
condition on its right end) is denoted as subdomain 2. 

We consider three different cases, and for the first two cases the trapezoidal parameter is taken 
as 7 = 0.1. For the chosen properties, a;™^^ = 1200 and 0^2'^^ = 1178.1, and the critical time steps 
for individual unconstrained subdomains are 2.083 x 10~^ s and 2.122 x 10"^ s, respectively. The 
upper bound for the stabilization parameter is 2.5 (see equation (I5.58P ). and the critical time steps 
for constrained subdomains are 

^^max ^ 2.083 X 10~^ (2.5 - 2a) , At^'"'' = 2.122 x 10~^ (2.5 - 2a) 

In the first case, we have taken a = 1, and the critical time steps for individual constrained 
subdomains will be 1.042 x 10~^ s and 1.061 x 10^^ s. We have taken the time step as At = 0.001 s, 
which satisfies the stability criteria (j5.58p . In Figure [19] we have shown the numerical results for 
temperature and rate of temperature using a = 1. The obtained numerical results are stable as 
predicted by the theory, and also there is no visible drift in the temperature or rate of temperature. 

In the second case, we illustrate that (as predicted by the theory) larger Baumgarte parameters 
affect the critical time step, and can make the algorithm unstable. To this end, we have taken 
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a = 2.6. We have again chosen At = 0.001 s as it satisfies critical time steps for individual 
unconstrained subdomains (and also for a = 0). But, according to the theory, the chosen time 
step may be unstable for a = 2.6. (Note that the energy method provides sufficient conditions 
for stability, and the bounds can be very conservative). In Figure [20] we have shown the obtained 
numerical results for temperature and its rate using the Baumgarte parameter a = 2.6. In Figure 

we have plotted the interface Lagrange multiplier against time for both a = 1 and a = 2.6. 
As one can see, the numerical results using a = 1 are stable whereas the results using a = 2.6 
are unstable. From the above discussion, we conclude that the Baumgarte parameter does affect 
the critical time step for 7 < 1/2, which is in accord with the theory. Also, from these numerical 
results one can see that the theoretical stability bounds obtained for the Baumgarte stabilized 
method using the energy method are rather sharp. 

In the third case, we have taken 7 = 1/2. For this case, the theory predicts that the Baumgarte 
stabilized domain decomposition method is unconditionally stable, and there is no restriction on 
the choice of the Baumgarte parameter a > 0. (Of course, to avoid ill-conditioned systems, which 
may give numerical results that are unreliable, one may not use a huge value for a.) We choose 
a = 2.6 (which is unstable under 7 = 0.1), and also a larger value of a = 10. As one can see from 
Figure [22l the numerical results are stable under 7 = 1/2 even for higher values of the Baumgarte 
parameter, which again agrees with the theory. 

7. CONCLUSION 

We have presented four variants of dual substructuring domain decomposition method for first- 
order transient systems. Two of them (the d-continuity and modified d-continuity methods) enforce 
the continuity of temperatures along the subdomain interface. The third one (the i>-continuity 
method) enforces the continuity of rate of temperature along the subdomain interface. The fourth 
method (the Baumgarte stabilized method) enforces the continuity of a linear combination of 
the temperature and rate of temperature along the subdomain interface. We have assessed their 
numerical stability and accuracy (in particular, drift in the constraints) under the generalized 
trapezoidal family of time integrators. 

Time stepping schemes from the trapezoidal family are primarily developed for solving ODEs, 
and the stability criteria widely used for solving first-order transient problems are (in general) not 
valid for transient domain decomposition methods (as the governing equations are DAEs instead of 
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ODEs). For example, under the d-continuity method (which enforces the subdomain equation and 
kinematic constraints at integer time levels) we have shown that all time stepping schemes with 
the trapezoidal parameter < 7 < 1/2 are unconditionally unstable. This is a surprising result, as 
in the case of ODEs the time stepping schemes with < 7 < 1/2 are stable (conditionally) and the 
mid-point rule with 7 = 1/2 is actually unconditionally stable. 

The modified d-continuity method enforces the subdomain equation at intermediate time levels, 
and kinematic constraints at integer time levels has been proposed. Under this domain decomposi- 
tion method, we have shown that all the quantities of interest (temperature, rate of temperature, 
and interface Lagrange multipliers) are all bounded provided that the time step is less than the 
critical time steps of individual unconstrained subdomains. 

The stability analysis of the i7-continuity method based on the energy method reveals that the 
rate of temperatures at integer time levels are bounded. However, solely based on the energy 
method, one cannot infer anything about the boundedness of the temperature and interface La- 
grange multiplier at integer time levels. One main drawback of the i^-continuity method is that one 
may get irrecoverable drift in the original constraint (i.e., continuity of temperature), which may 
make the numerical results unphysical. 

We have also performed stability analysis on the Baumgarte stabilization method for first-order 
systems, and derived simple bounds for the user-defined parameter that ensures stability. To our 
knowledge no prior work has derived bounds on the Baumgarte parameter for index-2 DAEs to 
ensure stability in a mathematically rigorous way. 

APPENDIX: IMPLEMENTATION DETAILS 

In this section we will outline a numerical algorithm for implementing the modified d-continuity 
method. One can easily modify this numerical implementation procedure for other domain decom- 
position methods presented in Section [H 

Using equations (|4.4p - (j4.5p . the governing equation for subdomain i can be written as 

(7.1) M^d^+^^ = {Mi - (1 - 7)AtK,) df) + At/S"+^^ + AtCjX^''+^^ 

where the matrix Mi := Mi + jAtKi. It is easy to see that the matrices Mi (i = 1, • • • ,S) are 
symmetric and positive definite. By substituting equation (17. ip into equation ()4.6p . the interface 
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problem can be written as 



(7.2) GA("+^) = R 



where the Schur complement operator G and interface vector R are defined as 



5 S 
j=l 1=1 



(7.3) G■.= Y,C^M,'CJ, R:=-Y,CiM,'}\ 



(7.4) := (M, - (1 - 7)AtK,) dS") + ft^^^ 

Remark 7.1. /t is easy to show that the Schur complement operator G is symmetric and positive 
definite (as, in this paper, we have assumed that there are no redundant constraints). Hence, the 
interface problem ()7.2p is solvable and has a unique solution (as G is invertible). 

Then a simple numerical algorithm for implementing the modified d-continuity method reads 

• Solve the interface problem (j7.2p to obtain the Lagrange multipliers A^"'^'^^ 

• Using the obtained interface Lagrange multipliers A'^'^''''''^ for each subdomain (i = 1, • • • ,5) 
obtain d-"^^^ by solving the subdomain equation (j7.ip . 

• Using equation solve for vf^^\ That is, v^^^'^^ = ^ (df^^'^ - 

• If desired, calculate the rates i^^"^^^ and Lagrange multipliers A*-""*"^^ at integer time levels 
using equations (j4.7p and (|4.8p . respectively. 

In the FETI method (see [9] , and also [TTj and references therein) , one employs an iterative solver (in 
this case, one can use the Preconditioned Conjugate Gradient method [11] as the Schur complement 
operator G is symmetric and positive definite) for the interface problem, and a direct solver for 
solving the governing equation for individual subdomains. 
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Figure 5. A typical 2D mesh used in the numerical convergence analysis. A compu- 
tational domain is divided into 4 subdomains of equal size. The mesh and subdomain 
interface are shown in the figure. 
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Figure 6 . Split degree of freedom system (7 = 1/4) : ^-continuity domain decom- 
position method, (a) v^^'^"'^ and v^^^^^ are bounded, (b) The quantities t;^^ and 

(n) 

are growing, which is in accord with the theory. 
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Figure 7. Split degree of freedom system (7 = 1/4): d-continuity domain decom- 
position method. The Lagrange multiplier A^"'"'"''') is bounded but A^"^ is growing, 
which is similar to the counterexample presented in Figure El 
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Figure 8 . Split degree of freedom system (7 = 1/4) : d-continuity domain decom- 
position method. The temperatures and are bounded for a while. But 



due to large (and growing) values of v 



in) In) 



and A^"); the whole numerical algo- 



(n) (n) 

rithm fails after sufficiently long time. Theoretically the quantities and are 
bounded Vn. But on a computer, both these values also blow up after some time 
because of finite precision arithmetic. At t = 0.69 s, the obtained numerical values 



are dj"^ « 0.1517, « -0.3483, and A^") « 8.6457 x 10^^. The exact values for 



temperature and Lagrange multipliers are, respectively, (P^^^ = (F^^'^^ ~ 0.0225 and 
^exact _ 0.1011. Also uote that in the earlier figures the values are plotted up to 
time t = 0.35 s, but in this figure, in order to see appreciable differences, we went 
up to i = 0.7 s. 
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Figure 9. Split degree of freedom system (7 = 1/4): Modified d-continuity method. 
In this figure, the rate of temperature (top) at integral time steps for both subdo- 
mains, and the interface Lagrange multiplier (bottom) at both integer and interme- 
diate time levels are plotted. As one can see, all these quantities are bounded, which 
agrees with the theory. Note that the rate of temperature and Lagrange multiplier 
at integer time steps are evaluated using equations (|4.7p and (|4.8p . respectively. 
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Figure 10. Split degree of freedom system (7 = 3/4): d-continuity domain decom- 
position method. The rate of temperature (top) at integer time levels are plotted 
for both subdomains (i.e., and "w^^)- The Lagrange multiplier (bottom) at both 
integer A(") and weighted A^"'"'"'^) time levels are also plotted. All these quantities 
are bounded, and the numerical results agree with the theory. 
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Figure 11. Split degree of freedom system (7 = 3/4): d-continuity domain decom- 

(n) (n) 

position method. For this case, both the temperatures d\ and are bounded. 
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Figure 12. d-continuity method: ID test problem. The trapezoidal parameter is 
taken as 7 = 1/4, and the time step is taken as At = 10~^ s. In this figure, the rate 
of temperature is plotted against the spatial coordinate at t = 0.001 s. As predicted 
by the theory, the d-continuity with 7 = 1/4 is not stable. 
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Figure 13. d-continuity method: ID test problem. The trapezoidal parameter is 
taken as 7 = 3/4, and the time step is taken as At = 0.001 s. The temperature 
(top) and rate of temperature (bottom) are plotted against the spatial coordinate 
at various time levels. 
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Figure 14. Modified d-continuity method: ID test problem. The trapezoidal pa- 
rameter is taken as 7 = 1/4 (top) and 7 = 3/4 (bottom), and the time step is 
taken as At = 0.001 s. The temperature is plotted against the spatial coordinate at 
various time levels. Under the modified d-continuity method, both the Baumgarte 
parameters 7=1/4 and 7 = 3/4, which is also observed in the numerical simulation 
of ID test problem. 
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Figure 15. d-continuity: 2D test problem. The trapezoidal parameter is taken as 
7 = 3/4, and the time step is taken as At = 0.001 s. The temperature profiles are 
shown for t = 0.05 s (top) and t = 0.5 s (bottom). 
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Figure 16. Modified <i-continuity: 2D test problem. The trapezoidal parameter is 
taken as 7 = 3/4, and the time step is taken as At = 0.001 s. The temperature 
profiles are shown for t = 0.05 s (top) and t = 0.5 s (bottom). 
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Figure 17. d-continuity method: Spatial numerical convergence analysis on ID 
and 2D test problems for temperature (top) and rate of temperature (bottom) using 
7 = 3/4. The analysis is performed at the time level t = 0.01 s. Standard linear 
and bilinear elements are used for ID and 2D problems, respectively. 
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Figure 18. Modified d-continuity method: Spatial numerical convergence analysis 
on ID and 2D test problems for temperature using 7 = 1/4 (top) and 7 = 3/4 
(bottom). The analysis is performed at the time level t = 0.01 s. Standard linear 
and bilinear elements are used for ID and 2D problems, respectively. 
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Figure 19. Baumgarte stabilized DD method. The Baumgarte parameter is taken 
as a = 1, trapezoidal parameter as 7 = 0.1, and the time step as 0.001 s. Each 
subdomain is divided into 10 equal elements. The obtained numerical solutions for 
temperature (top) and rate of temperature (bottom) are plotted against x at time 
levels 0.25 s to 2.00 s with time increments of 0.25 s. 
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Figure 20. Baumgarte stabilized DD method. The Baumgarte parameter is taken 
as a = 2.6, trapezoidal parameter as 7 = 0.1, and the time step as 0.001 s. Each 
subdomain is divided into 10 equal elements. The obtained numerical solutions at 
t = 0.8 s for temperature (top) and rate of temperature (bottom) are plotted against 

X. 
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Figure 21. Baumgarte stabilized DD method. The trapezoidal parameter as 7 = 
0.1, the time step as 0.001 s, and each subdomain is divided into 10 equal elements. 
In this figure we have plotted the interface Lagrange multiplier against time for the 
Baumgarte parameters a = 1 (top) and a = 2.6 (bottom). 
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Figure 22. Baumgarte stabilized DD method. The trapezoidal parameter is taken 
as 7 = 1/2, and the time step as 0.001 s. Each subdomain is divided into 10 equal 
elements. The temperature profiles for a = 2.6 (top) and a = 10 (bottom) are 
plotted against x at time levels 0.25 s to 2.00 s with time increments of 0.25 s. 
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